Ecogenomics reveals viral communities across the Challenger Deep oceanic trench

Despite the environmental challenges and nutrient scarcity, the geographically isolated Challenger Deep in Mariana trench, is considered a dynamic hotspot of microbial activity. Hadal viruses are the least explored microorganisms in Challenger Deep, while their taxonomic and functional diversity and ecological impact on deep-sea biogeochemistry are poorly described. Here, we collect 13 sediment cores from slope and bottom-axis sites across the Challenger Deep (down to ~11 kilometers depth), and identify 1,628 previously undescribed viral operational taxonomic units at species level. Community-wide analyses reveals 1,299 viral genera and distinct viral diversity across the trench, which is significantly higher at the bottom-axis vs. slope sites of the trench. 77% of these viral genera have not been previously identified in soils, deep-sea sediments and other oceanic settings. Key prokaryotes involved in hadal carbon and nitrogen cycling are predicted to be potential hosts infected by these viruses. The detected putative auxiliary metabolic genes suggest that viruses at Challenger Deep could modulate the carbohydrate and sulfur metabolisms of their potential hosts, and stabilize host’s cell membranes under extreme hydrostatic pressures. Our results shed light on hadal viral metabolic capabilities, contribute to understanding deep sea ecology and on functional adaptions of hadal viruses for future research.


R1:
We thank this reviewer for the overall positive comments, and we appreciate his/her concern. However, we need to that point out that the viral communities reported by Zhao et al., (2022) are from the upper slope of the trench, and specifically, from surficial sediments (0-18 cm) at four sampling sites between 5.4 to 6.7 km water depth (see Table S1 by Zhao et al., 2022). In contrast, our study reports data on viral communities from 13 sampling sites that cover both slopes and the bottom-axis of the Challenger Deep, targeting primarily hadal depths below 6 km, with the deepest sampling depth at ~11 km. Although Challenger Deep is still understudied, we know from literature that the funneling effect and geographic isolation creates heterogeneity between sites along the slopes of the trench and its bottom axis, shaping prokaryotic microbial communities and diversity (e.g., Luo et al., 2017Zhou et al., 2022. Because of this we believe that even intra-trench comparisons of microbial and viral communities (e.g., slope vs. slope sites or bottom-axis vs. bottom-axis sites) will produce new findings. The inclusion of a greater number of sites as well as sites that extend down to ~11km in our study therefore does not just extend the findings of viruses in the Zhao et al. 2022 study. Nonetheless, to address this reviewer's concern, we compared our viral data to the data reported by Zhao et al. 2022, and found that 98% are new species (<95% identity in 85% of sequence length) and 76% are new genera (estimated by vcontact2). We amended one sentence in our abstract to say "Here, we collected 13 sediment cores from slope and bottom axis sites across the Challenger Deep (down to ~11 km depth), and identified 1,628 previously undescribed viral operational taxonomic units at the species level." And in the main text in lines 159-163 we now say: "Our CD vOTUs were also novel when compared with viromes identified at the hadal slope sediments of the Mariana Trench 48 . Specifically, 98% of our CD viral contigs were new species (< 95% identity in 85% of sequence length) and 76% of them were new genera (estimated by vcontact2), when compared with the identified viruses from the upper slope (5.4-6.7 km depth) of the trench 48 ." For reviewer's convenience the papers referred to above are here: Luo M, Gieskes J, Chen L, et al. Provenances, distribution, and accumulation of organic matter in the southern Mariana Trench rim and slope: Implication for carbon cycle and burial in hadal trenches [J]. Marine Geology, 2017, 386: 98-106. Zhou Y, Mara P, Cui G, et al. Microbiomes in the Challenger Deep slope and bottom-axis sediments [J]. Nature Communications, 2022, 13(1): 1-13. Zhao J, Jing H, Wang Z, et al. Novel Viral Communities Potentially Assisting in Carbon, Nitrogen, and Sulfur Metabolism in the Upper Slope Sediments of Mariana Trench. mSystems. 2022;7(1):e0135821. Specific comments: Q2: As you use the viral contigs after manual curation (1628 vOTUs), I'd like to know the performance comparison of the tools about the high confidence viral contigs. Was there any software identifying all high confidence vOTUs?

R2:
We performed laborious manual curation of the viral contigs to avoid false-positive predictions, and PPR-meta analysis to identify the most confidently-assigned viral contigs by the twelve tools incorporated into our pipeline (Figure 2a). We now provide Supplementary Table 4 that summarizes the screening criteria used for identifying the viral contigs, and the cut-offs used for each tool. We didn't a single individual software to identify the high confidence vOTUs. We now include the term "putative" viral contigs in various parts of the paper after reviewer's 2 request.

Q3:
Line 37 slope-axis should be mentioned in abstract, otherwise we do not know what the bottomaxis was compared with.

R3:
This reviewer is right. We rephrased accordingly and now reads "Community-wide analyses revealed distinct viral diversity across the trench which is significantly higher at the bottom-axis, when compared to slope sites. In silico predictions indicate viral infection of key prokaryotes involved in hadal carbon and nitrogen cycling". Please see lines 36-38.

Q4:
Line 83-113 Please briefly introduce what you did in this paragraph, not an extended edition of abstract.

R4:
Lines 82-113 describe environmental features of Challenger Deep (e.g., geochemistry, in situ temperature/pressure conditions, carbon distribution) and summarize what is known for hadal microbial communities. This information is not described/discussed in the abstract and we consider it of importance to the reader. We also believe that lines 82-84 and 105-113 that describe the sampling sites and the aims of this study do not overlap with what is already discussed in abstract.

Q5:
Line 120 I'd like to know the relative transcription level of viruses compared to hosts (something like certain prokaryotic housekeeping gene?), even though the percentage might be very low.
R5: This is a really interesting comment. Our data allowed us to predict hosts for a small fraction of the CD viral community (14/1628 vOTUs; see lines 230-232), and we also have only three metatranscriptome libraries (6-9, 12-15, 18-21 cmbsf) from one sediment core collected at ~11 km depth. Most genes were not mapped to by metatranscriptome reads in our libraries. For those transcripts that mapped, a significant fraction was characterized as genes of unknown function (16% of the total mRNA pools; please see Supplementary Data 4, in Zhou et al., 2022). We do not consider this surprising since the available mRNA data from hadal trenches are limited, and we also know that RNA can be readily degraded especially when extensive recovery times of samples occur (~5 hrs from 11 km depth to surface; this study). This would make it challenging even for genes from the core genome (e.g., housekeeping genes) to be reliably used for any inter-sample comparison. To address this comment, we listed the number of host genes that were mapped by metatranscriptome reads in the

R7:
We agree with the reviewer on this. However, providing percentages might not be as precise or widely accepted for estimating viral activity especially when bulk metatranscriptomic pools are used.
Nonetheless, to follow this reviewer's advice, we now provide Supplementary Table 6 that shows the different % of genes mapped by each vOTU using metatranscriptome reads, and the different percentages of potentially active genes in the viral contigs (see below). We also rephrased the text (see lines 143-145) to clarify that we use 20% as threshold for identifying potential active viral contigs. Supplementary Q8: Line 227 Please include the relative abundance of the potential host MAGs. The discussion about the influence of lytic viruses in hadal carbon and nitrogen cycling are based on the widespread and high abundance of these potential hosts.

R8:
We added a sheet to Supplementary Table 7 to show the relative abundance of the potential host MAGs in each metagenome (please see sheet "Relative abundance"), and we also provide this information in the text. Now on line 235-239 it reads: "These taxa include heterotrophs (e.g., Proteobacteria) and chemoautotrophs (e.g., Thaumarchaeota, Planctomycetota) involved in nitrogen and carbon cycling whose taxonomic signatures were abundant in CD sediments 39 , but with different relative abundances (7% to 43%) across the discrete sampling sites (bottom-axis vs. slope) (Supplementary Table 9)." See lines 572-576 for the method for estimate of the relative abundance of host MAGs.

Q9:
Line 277 This conclusion is too strange considering you only predicted the hosts of 14 vOTUs. R9: We agree with the reviewer and we rephrased accordingly. Now on line 280-283 it reads: "The predicted potential prokaryotic hosts for the 14 vOTUs may suggest that CD viruses target specific prokaryotic hosts in these CD sediments, however, this requires further investigation considering that our host predictions were accomplished for ~1% of the viral population that we identified.
Q10: Line 456 Annotated using which database? eggnog? R10: Yes, we used the eggNOG database for annotation. We added text to clarify this and it now reads: " 2. We annotated the putative virus contigs using the eggNOG database." See lines 463-464.

Q11:
Line 471 Please include the standard considering a viral contig as contamination. The coverage or the covered percentage?

R11:
We added text to clarify. Now it reads (lines 482-484): "Viral contigs mapped with ≥ 1 reads from a blank control were considered potential contaminants. This resulted in the removal of ten viral contigs that were excluded from further analysis." Q12: Line 477 Considering you performed vConTACT2, please include the taxonomy assignment results of this software.
R12: vConTACT2 clustered/classified only one viral contig using its reference database. This did not change the results of the taxonomy annotation that we already had.

Q13:
Line 553 "a viral contig database"? R13: We apologize for this vague statement. We were referring to all CD viral contigs. We rephrased lines 578-580 as follows: "To calculate the coverage (sequencing depth) of each viral contig, clean and qualified reads from each sample were mapped against all CD viral contigs using BWA (v 0.7.17) and sorted with samtools (v1.9).
Q14: Line 558 The threshold value of 10% is too low, which might cause too many false data. Please consider the threshold value of vOTUs clustering (85%).

R14:
We thank the reviewer for this suggestion. However, the sequencing depths of CD viral contigs is not high, and for this reason we would like to be as stringent as possible, and avoid false-positive outcomes. Our viral contigs could only recruit every few viral reads from the bulk metagenomes.
Using higher thresholds, as suggested by the reviewer, would calculate most viral contigs as zero coverages creating false-negative results due to the low sequencing depths. Below you can find three figures that we generated for this reviewer, using the threshold that he/she suggested (85%) (Fig. a), a 50% threshold (Fig. b) and the 10% threshold (Fig. c) which is used in this study. Considering that almost all (1,622/1,628) of our viral contigs are > 10kb, a threshold of 10%, will result in many/most viral contig to still have > 1kb region covered by reads (identity > 95%). This reduces false-positive outcomes. Roux et al. 2017, reports that increasing thresholds (read mapping identity percentage and length of contig covered), progressively decreases the sensitivity of the analysis, and the false discovery rate (defined as percentage of contigs recovered, that were not part of the initial community). This is something that we would like to avoid with our CD data. However, we report both alpha and beta diversities (plus other indices Chao, Shannon etc), so even with low coverage (viral contigs > 10% of their length covered by metagenome reads) we can still characterize fairly well the viral diversity of CD. We know that alpha diversity can be highly variable when samples are significantly under-sequenced, but beta diversity trends can be recovered even when sequencing depth is highly variable (

Q15:
We added commas where necessary, following reviewer's suggestion.

Q1:
The authors took on the endeavor to detect, characterize and compare viruses recovered from deep ocean trenches. I appreciated all the efforts to study these under-sampled environments that will advance our understanding of environmental viruses and their ecological functions. Overall, the topic is important, and the results are exciting to read. However, I do have some major concerns. 1) More thoughts are needed in methodology: example 1: the authors applied one unpublished or not peer-reviewed workflow to detect viral contigs. Although the tools mentioned in the workflow are widely used, the cut-offs and ways of sorting the results are encrypted in the workflow; example 2: suspicious methods to screen putative viral contigs such as 'at least 70% of proteins in the contig were assigned as 'hypothetical protein', 'unknown function'; example 3: identify lytic viruses using VIBRANT that can lead to misinterpretation of the results.
2) The authors need to be careful when citing references to support your discussion. Examples: citing soil viruses for supporting low lysogeny; citing thawed permafrost papers for permafrost ecosystem; citing Paez-Espino et al. (2017) for the suspicious method of screening putative viral contigs.

R1:
We thank for the reviewer for finding our results exciting to read. We address this reviewer's concerns on our point-by-point responses.
Q3: Line 57-61 and line 69-73 please consider separating the long sentence into two. R3: We rephrased lines 57-61 (now new lines 57-61) following reviewer's suggestion. Now it reads: "Viruses show high abundances in marine sediments (10 7 -10 10 particles g -1 of dry sediment). Yet, viral particles bind firmly to sediments due to electrostatic, "van der Waals" and hydrophobic interactions, which complicate their separation and enumeration from the surrounding sediment matrix 21" . We did the same for lines 69-73 (now new lines 69-71). Now it reads: "Among prokaryotes, Thaumarchaeota and other archaeal lineages in deep-sea sediments, are reported to be more susceptible to viral infections compared to bacterial taxa 27 ". We have deleted lines 71-72 (see Q4) to avoid redundancy. The notion in lines 71-72 was already addressed in lines 67-69: "The viral shunt in abyssal and hadal realms is estimated to contribute 35% of labile carbon in those habitats and is believed to sustain the sediment microbiota in hadal sediments by providing easily degradable carbon" Q4: Line 71-72, 'the fast decomposition of released viruses following prokaryotic infections'. Do you mean 'the fast decomposition of viruses released by the lysed prokaryotic host cells'? if so, 'prokaryotic infection' is misleading. Please consider re-writing it for clarity.
Q5: I would encourage the authors to enclose the sequencing and assembly statistics (e.g., qualityfiltered reads, Numbers of contigs/scaffolds, N50, numbers of reads that contribute to viral contigs etc.) in the supplementary file for studies using metagenomes/metatranscriptomes.

R5:
Thanks for the suggestion. Please find new Supplementary Table 3, that now includes the sequencing and assembly statistics.

Q6:
Line 121, what is 'T3L11'? why this was selected? If there is no particular reason, please consider editing it into 'from one of the sediment cores (T3L11, 10,908 m)'.

R6
: T3L11 was the deepest site sampled during our cruise in Challenger Deep. Also, T3L11 is the only sediment core for which we have metatranscriptome data. We corrected the text following the reviewer's suggestion and now it reads: "We also generated three metatranscriptome libraries from one of the bottom-axis sediment cores (T3L11: 10,908 m; 6-9, 12-15, 18-21 cmbsf) to gain insights into potential viral activities". Please see lines 120-122.

Q7:
It is a bit worrying that the authors used viral detection workflow that is not published or peer reviewed (Marquet, Mike, et al. "What the Phage: A scalable workflow for the identification and analysis of phage sequences." bioRxiv (2020).).

R7:
We understand this reviewer's concern. The pipeline by Marquet et al., combines 12 tools for phage annotation and identification. We used it primarily to check if viral contigs were present in our data set, and if yes, if they had sizes of > 10kb. Considering many/most viral contigs in CD are novel, it would be beneficial to know the performance of different viral prediction tools for these data. However, because this is not a peer-reviewed pipeline, and also had a highly variable prediction quality, we were also skeptical and wanted to be as stringent as possible. We therefore performed rigorous manual curation of the data as we explain in lines 129-136. This removed more than 80% of the generated viral contigs from further downstream analysis. The predicted viral contigs used for this study are ~16% (1628/9889) of those initially identified by the pipeline and are retained using consensus metrics of > 95% identity and > 85% coverage. Following this reviewers' suggestion on Q20, we also now include Supplementary Table 4 that summarizes the screening criteria for identifying the viral contigs, the cut-offs used for each tool and we include the term "putative" viral contigs to be more careful how we refer to them in our study.
We believe that what we report are reliable and uncontaminated data, and that in our effort to be as cautious as possible we have potentially excluded many real viral contigs in the course of our analyses.
Q8: Line 131-134, a total of 1628 contigs were detected and they were clustered in 1628 vOTUs with sequences longer than 10 kb and another six vOTUs with sequences shorter 10 kb? These sentences are confusing.

R8:
We apologize for that. We rephrased to avoid any misunderstanding. Now it reads: "Overall, 1,622/1,628 vOTUs were > 10 kb while six had sizes less than 10 kb ( Q9: Line 135, 'due to the removal of host regions from the proviruses' R9: Please see our R8 response. We only identify viruses from contigs > 10kb. After viral identifications, we used checkv to estimate viral completeness, which will also remove host regions from proviruses.

Q10:
Line 137, what does it mean, 'at least for vOTUs with closely related reference genomes'? if certain criteria were applied to assess the vOTUs using checkV, please write in a full sentence. If not, please consider removing it for clarity.

R10:
We apologize for this vague statement. We rephrased accordingly and now it reads: "The degree of completeness and contamination of the CD vOTUs was estimated by comparing the sequences using CheckV 44 against a large database of environmentally-diverse, and complete viral genomes.". See lines 136-138.
Q12: Line 164, 'estimated relative abundances' ? because the reads coverage was calculated by the reads that were mapped to viral contigs relative to the ones that were not. Reads coverage is just an estimate.

R12:
We agree with this reviewer and we corrected this as suggested. See line 168: "The estimated abundances of vOTUs" Q13: Line 169, do we know the taxonomy of 'T1L10_NODE_10823'? any close relatives in the reference database?
R13: This would be ideal but unfortunately, we do not know the taxonomy of T1L10_NODE_10823. Also, T1L10_NODE_10823 had no close relatives in the reference databases, and was not affiliated with the identified Thaumarchaeota viruses.
Q14: Line 180, it is worrying to classify the viral contigs that were not detected as 'prophages' or 'potential temperate viruses' are lytic viruses, although the authors acknowledged the risks of overestimation.

R14:
We understand this concern and for this reason we rephrased all text between lines 187-198 and we replotted Figure 2b (now Figure 2e). This figure now describes viral contigs as "Lysogenic" and "Unassigned" (see below). Also, the rephrased text now includes the suggestions of Q15-Q16 by this reviewer. Please see lines (187)(188)(189)(190)(191)(192)(193)(194)(195)(196)(197)(198): "Our results indicated that 1,541 viral contigs (95%) in CD viromes were not assigned to either a lytic or lysogenic lifestyle (Fig. 2e, "undetermined"). It is possible that many/most of these "undetermined" viral contigs belong to viruses that have a lytic lifestyle in hadal depths. This would be consistent with studies of viral communities from surficial sediments collected in different deep-sea oceanic settings (Arctic, Atlantic, Pacific Oceans and Mediterranean Sea; > 1,000 m water depth) that report high viral lysis rates 27 . With regard to lysogeny, it was predicted only in 5% of the CD sediment viromes. This differs from deep-sea sediments that showed lysogeny as a more common potential viral lifestyle (e.g. Baltic Sea; ~19% on average) 25 but is more in line with the prediction results that we obtained for deep-sea cold seep sediments (7%) 15 and ocean seawater viromes (3%) 7 using VIBRANT 49 . Nonetheless, our arguments need to be interpreted with caution considering that 95% of viral contigs were not assigned as lytic or lysogenic. Now, we only predict lysogenic viruses and we re-made figure 2e to reflect this change. All other viruses are classified as "undetermined."

Fig 2e
Q15: Line 189, please double check the two soil references to support the 5% of lysogeny rate in soil. In soil, the rate is relatively more accurately estimated using induction assay. Lysogeny is quite prevalent in soil. Soil may be not a comparable environment in this case. (ref: Incidence of lysogeny within temperate and extreme soil environments; Prevalence of Lysogeny among Soil Bacteria and Presence of 16S rRNA and trzN Genes in Viral-Community DNA).

R15:
We appreciate this important comment from this reviewer. We have now removed soil samples and focused our comparisons only on viral data reported from oceanic settings (e.g., water column, deep-sea cold seep sediments) Please see our R14 for revised text.
Q16: Line 189, could you show me the content in reference 7 indicate the lysogeny percentage of 3%?

R16:
We apologize for this misunderstanding. Reference 7 cites the data (and not the 3%) that we used to perform the same prediction analysis of viral lifestyle that we performed for our vOTUs, and we compare the findings. Please see our R14 for revised text.
Q17: Line 210, please use 'viral' instead, as 'virome' refers to the viral fraction that are experimentally enriched from the environmental samples.
Q18: Line 255-263, I would suggest shortening this discussion since there are lots of speculations based on unconfident analysis of lytic viruses.

R18:
We agree with this reviewer and we shortened this part. The new text on lines 261-266 now reads: "Based on our analyses, lysogeny is a less likely lifestyle (5% assigned) in our identified CD viral contigs. Yet, the inability to assign lifestyle to the majority of the viral contigs (95%) might underestimate the importance of lysogeny, while at the same time preventing us from predicting the lytic viruses in CD. We suggest that lytic infections (if occurring) might be important and affect available nutrient pools across the V-shaped Challenger Deep (bottom-axis vs. slopes sites)." Q19: Line 321-324, how the modeled protein structure can predict the potential enzymatic activity? How the structure is similar to the reference with validated activity? can some of the residues be aligned to the active sites of these reference structures if any?
R19: All these are great questions. We have now added text in the Methods regarding the modelling of the protein structure that we report on our discussion of putative AMGs (lines 536-543, see below).
We also clarify in the discussion how the structure predictions can provide information on the potential enzymatic role that we believe might be useful for the reader (lines 328-339 and below).
Overall, in our study we used the web-based Phyre2 services for protein structure prediction Structural homologies are estimated against publicly available proteins whose function is known and structure has been solved using crystallography or other appropriates techniques. We could have performed additional analysis e.g., PROSITE-Expasy, to describe in more depth the active domains and functional sites at the residue level in our predicted CysC proteins, however we believe that the results we report with Phyre2 are sufficient for the purpose of this study.
Below we include Figure 1  For this reviewer our added discussion on lines 328-339: "The distinct phylogenetic results and the moderate similarity of the CD Cys proteins to those that are publicly available, prompted us to perform protein structure prediction for the CysC protein from the viral contig T1B8_NODE_1222 (Fig. 6c). We used the web-based Phyre2 tool that predicts protein structure and function using homology with known proteins available in protein data banks 59 (see Methods Also, the added text in the Methods on lines 536-543 now reads: "The structure prediction for CysC protein was performed with the web-based Phyre2 tool 59 . Structural homologies were analyzed using models generated by Phyre2 using a confidence threshold of > 98%, and identity threshold of > 29%. The accuracy of the models constructed using Phyre2 is described as extremely high when the sequence identity is above 30-40%. However, lower sequence identities can be equally accurate and useful as long as the confidence threshold is high, which was the case in our examined CysC proteins. The functional domain for CysC was identified and annotated by SMART 113 . Figure 1 (from Kelly et al., 2015).

R21:
Thank you for this question. We will elaborate to avoid potential misunderstandings. Lines 458-481 explain the criteria we used for screening viral contigs; however we rephrased this part for clarity (see lines 463-467). Overall, our criteria for identifying the viral contigs were: 1) contig size (>10kb), 2) presence of 2 or more hallmark viral genes in the contig and 3) absence of any prokaryotic signature in the contig; in case the contig contained prokaryotic-specific genes it was removed from further analysis. We don't use their pipeline (Paez-Espino et al., 2017), and we apologize if this was interpreted this way. We cited Paez-Espino et al., 2017 because we wanted to show that the authors retained putative viral contigs that mostly contained genes of unknown and hypothetical function as long as 1) contigs included hallmark viral genes and 2) had absence of plasmid or microbial signature gene sequences. Also, Paez-Espino et al., (2017) screened all DNA metagenomic contigs that were more than 5 kb in size, while the authors use as a filter (Filter 1 out of 3) the following: "metagenomic contigs that had at least 5 hits to viral protein families; AND Total number of genes covered with KO terms on the contig ≤20%. AND Total number of genes covered with Pfams ≤40%; AND Total number of genes covered with viral protein families ≥10%."). These criteria identified contigs as viral contigs in their study, although they contained high percentage of unknown function genes. Also, we don't find it surprising that contigs contain genes encoding hypothetical proteins or proteins of unknown function since publicly available databases are not enriched with viral data, especially coming from deep-sea and hadal trench habitats. Further, culture-based experiments that could be more precise in assigning functions to proteins exist primarily for fungi and bacteria from Challenger Deep (at least to our knowledge). Besides Paez-Espino et al (2017), we also now cite Gao et al., 2020 (Gao, SM., Schippers, A., Chen, N. et al. Depth-related variability in viral communities in highly stratified sulfidic mine tailings. Microbiome 8, 89 (2020); https://doi.org/10.1186/s40168-020-00848-3), that also retained viral contigs that had a total number of genes assigned as "unknown" (annotated with eggNOG v5.0.0 database) for ≥ 80% of the total number of genes on viral scaffolds (> 10kbp in size) or scaffolds were enriched in hypothetical proteins.

Q22:
Missing the method section of mapping metatranscriptome to viral contigs.

Figure 6b
Maximum-likelihood phylogenetic tree of CD CysC proteins. The CysC proteins predicted in CD viral genomes were used to construct a phylogenetic tree using homologous CysC proteins deposited in the eggNOG database (V5.0) and publicly available viromes. CysH proteins were used as an outgroup. We also included two CysC homologs from the Uniport database (in blue) with experimental evidence of function at the protein level. Bootstrap values (1,000 replicates) ≥ 70% are indicated at nodes.